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[i]  How  to  reduce  the  horizontal  pressure  gradient  error  is  a  key  issue  in  terrain-following 
coastal  models.  The  horizontal  pressure  gradient  splits  into  two  parts,  and  incomplete 
cancellation  of  the  truncation  errors  of  those  parts  cause  the  error.  Use  of  the  finite  volume 
discretization  leads  to  a  conserved  scheme  for  pressure  gradient  computation  that  has 
better  truncation  properties  with  high  accuracy.  The  analytical  coastal  topography  and 
seamount  test  cases  are  used  to  evaluate  the  new  scheme.  The  accuracy  of  the  new  scheme 
is  comparable  to  the  sixth-order  combined  compact  scheme  (with  an  error  reduction  by  a 
factor  of  70  comparing  to  the  second-order  scheme)  with  mild  topography  and  much 
better  than  the  sixth-order  combined  compact  scheme  with  steep  topography.  The 
computational  efficiency  of  the  new  scheme  is  comparable  to  the  second-order  difference 
scheme.  The  two  characteristics,  high  accuracy  and  computational  efficiency,  make 
this  scheme  useful  for  the  sigma  coordinate  ocean  models.  Index  terms:  4255 
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1.  Introduction 

[2]  In  regional  oceanic  (or  atmospheric)  prediction  mod¬ 
els,  the  effects  of  bottom  topography  must  be  taken  into 
account  and  usually  the  terrain-following  sigma  coordinates 
should  be  used  to  imply  the  continuous  topography.  In 
sigma  coordinates  the  water  column  is  divided  into  the 
same  number  of  grid  cells  regardless  of  the  depth.  Consider 
2D  problems  for  mathematical  simplification.  Let  (x,  z)  be 
the  Cartesian  coordinates  and  (x,  a)  be  the  sigma  coordi¬ 
nates.  The  conventional  relationships  between  z  and  sigma 
coordinates  are  given  by 


where  t]  is  the  surface  elevation.  Both  z  and  a  increase 
vertically  upward  such  that  z  =  r|,  a  =  0  at  the  surface  and 
=  —1,  z  =  —H  at  the  bottom.  The  horizontal  pressure 
gradient  becomes  difference  between  two  large  terms 

?P  =  ?P _ 

dx  dx  H  +  a\  dx  dx)  do'  {) 

that  may  cause  large  truncation  error  at  steep  topography 
[e.g.,  Gary,  1973;  Haney,  1991;  Mellor  et  al.,  1994; 
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McCalpin,  1994;  Chu  and  Fan,  1997,  1998,  1999,  2000, 
2001;  Song,  1998], 

[3]  Several  methods  have  been  suggested  to  reduce  the 
truncation  errors  to  acceptable  levels:  (1)  smoothing  topog¬ 
raphy  [Chu  and  Fan,  2001],  (2)  subtracting  a  mean  vertical 
density  profile  before  calculating  the  gradient  [Gary,  1973; 
Mellor  et  al.,  1994],  (3)  bringing  certain  symmetries  of  the 
continuous  forms  into  the  discrete  level  to  ensure  cancella¬ 
tions  of  these  terms  such  as  the  density  Jacobian  scheme 
[e.g.,  Mellor  et  al.,  1998;  Song,  1998;  Song  and  Wright, 
1998],  (4)  increasing  numerical  accuracy  [e.g.,  McCalpin, 
1994;  Chu  and  Fan,  1997,  1998,  1999,  2000,  2001],  (5) 
changing  the  grid  from  a  sigma  grid  to  a  z  level  grid  before 
calculating  the  horizontal  pressure  gradient  [e.g.,  Stelling 
and  van  Kester,  1994],  Kliem  and  Pietrzak  [1999]  claimed 
that  the  z  level  based  pressure  gradient  calculation  is  the 
most  simple  and  effective  means  to  reduce  the  pressure 
gradient  errors.  After  comparing  to  other  schemes,  Ezer  et 
al.  [2002]  show  the  favorable  performance  of  the  latest 
polynomial  schemes.  Recently,  Shchepetkin  and  McWil¬ 
liams  [2002]  design  a  pressure  gradient  algorithm  with 
splines  that  achieves  more  accurate  hydrostatic  balance 
between  the  two  components  and  that  does  not  lose  as 
much  accuracy  with  nonuniform  vertical  grids  at  relatively 
coarse  resolution. 

[4]  Using  the  finite  volume  integration  approach  [Lin, 
1997],  an  extra  hydrostatic  correction  term  is  added  to  the 
ordinary  second-order  scheme  for  reducing  the  horizontal 
pressure  gradient  error.  With  the  extra  term  the  discretiza- 
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Figure  1.  Finite  volume  discretization  with  staggered  grid 
in  a  terrain-following  coordinate  system. 


tion  scheme  is  called  the  hydrostatic  correction  (HC) 
scheme.  In  this  study,  we  describe  the  physical  and  math¬ 
ematical  bases  of  the  HC  scheme  and  its  verification.  The 
outline  of  this  paper  is  as  follows:  Description  of  the 
horizontal  gradient  in  finite  volume  is  given  in  section  2. 
The  second-order  scheme  is  given  in  section  3.  Hermit 
polynomial  integration  and  hydrostatic  correction  are 
depicted  in  sections  4  and  5.  Error  estimation  and  seamount 
test  case  are  given  in  sections  6  and  7.  In  section  8,  the 
conclusions  are  presented. 

2.  Horizontal  Pressure  Gradient  in  Finite  Volume 

[5]  Let  the  flow  field  change  in  x-z  plane  only  (Figure  1). 
A  finite  volume  (trapezoidal  cylinder)  is  considered  with  the 
length  of  Ly  (in  the  y  direction)  and  the  cross  section 
represented  by  the  shaded  region  (trapezoid)  in  Figure  1. 
The  resultant  pressure  force  (F)  acting  on  the  finite  volume 
is  computed  as  follows: 

F  =  Ly  j)  pads  (3) 

c 

where  p  is  the  pressure,  C  represents  the  four  boundaries;  n 
denotes  the  normal  unit  vector  pointing  inward;  and  ds  is  an 
element  of  the  boundary.  The  contour  integral  is  taken 
counterclockwise  along  the  peripheral  of  the  volume 
element.  The  pressure  force  exerts  on  the  four  boundaries 
of  the  finite  volume  with  pw  pe,  pu,  and  pt  on  the  west,  east, 
upper,  and  lower  sides.  The  horizontal  (Fx)  and  vertical  (Fz) 
components  of  the  resultant  pressure  force  are  computed  by 


where  points  1,  2,  3,  and  4  are  the  four  vertices  of  the  finite 
volume.  The  hydrostatic  balance  is  given  by 


Fz  -  gAm,  (6) 

where  g  is  the  gravitational  acceleration,  Am  is  the  mass  of 
the  finite  volume.  Equation  (6)  states  that  the  vertical 
component  of  the  resultant  pressure  force  acting  on  the 
finite  volume  exactly  balances  the  total  weight  of  the  finite 
volume. 

[6]  For  a  Boussinesq,  hydrostatic  ocean  model,  the  pres¬ 
sure  field  is  calculated  by 

0 

P=Pam  +  PoCT  +  g  j  P(*>z' -  »  (?) 

z 

where  patm  is  the  atmospheric  pressure  at  the  ocean  surface, 
p0  is  the  characteristic  density,  and  tj  is  the  surface  elevation. 
Substitution  of  (7)  into  (5)  leads  to 


2  0 


4  0 


Fz  =  gLy  (//  p(x,z!  ,t)dz!  dx  +  II  p(x,z!  ,t)dz'dxj 
=  gLy  J  J  p(x,  z' ,  t)dJdx  =  gAm, 


(8) 


where  AS  is  the  area  of  the  trapezoid  (Figure  1)  computed 
by 

AS  =  (xi+i  -  xfifaji  +  zt+ijc  -Zijc+i  -Zi+\j[+i)tZij[=Hrak. 

(9) 

Equation  (8)  indicates  that  the  finite  volume  discretization 
guarantees  the  hydrostatic  balance  in  Boussinesq,  hydro¬ 
static  ocean  models.  Using  equation  (4)  the  horizontal 
pressure  gradient  is  computed  by 


dp  _  Fx 
dx  ~  LyAS 


which  has  the  flux  form  (i.e.,  the  conserved  scheme).  In  the 
conserved  scheme,  the  pressure  gradient  for  any  grid  cell  is 
computed  from  the  summation  of  pressure  exerted  on  the 
four  sides  of  the  cell.  For  the  whole  domain  integration,  the 
pressure  at  any  grid  cell  side  not  at  the  domain  boundary  is 
used  twice  with  opposite  sign  (canceling  each  other).  Thus 
equation  (10)  is  called  the  conservation  scheme.  Finite 
difference  schemes  in  the  sigma  coordinate  such  as  second- 
order  central  difference  scheme  popularly  used  in  ocean 
models  as  well  as  the  recently  developed  spline-based 
scheme  [Shchepetkin  and  McWilliams ,  2002].  Since  design 
of  conservation  scheme  is  a  key  issue  in  numerical  modeling, 
the  discretization  (equation  (10))  may  improve  the  pressure 
gradient  computation  in  terrain-following  ocean  models. 


3.  Second-Order  Scheme 


[7]  In  ocean  models  with  the  staggered  grid  (Figure  1), 
velocity  is  evaluated  at  the  center  of  the  volume  and 
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pressure  is  put  at  the  four  vertices.  With  given  values  at 
vertices,  we  have  several  methods  to  compute  four  integra¬ 
tions  in  the  right-hand  side  of  (equation  (10)).  The  mean 
value  theorem  leads  to, 


L  -  5 

J  Ptdz  =  Pi  (Z/+ 1  Jc+l  Z/ji+l ) ,  J  pedz  =  pe  (z,+ 1  .k  —  Z|'+ 1  .k  •  1 ) , 

1  2 
4  1 

J Pidz  —  pu(ziji  —  J  pwdz  =  pw  (zifi+i  —  Zijt) ,  (11) 


where  pi,pe,pu,pw  are  the  mean  values  of  pressure  at  the 
four  sides  of  the  trapezoid.  The  horizontal  pressure  gradient 
with  the  finite  volume  consideration  is  given  by 

~Kx  =  AS  [P'(z,+1’*+1  ~ ZA+0  +Pe{Z‘+l*  ~  Z'+l,*+l) 

+Pu(z>,k  -  zi+ 1,*)  +pw(z,Mi  -  zitk)] ,  (12) 


For  the  second-order  staggered  grid ,  pi,pe,pu,Pw,  are  taken 
as  the  arithmetic  means  of  pressure  at  the  two  vertices, 


_ ,  Pi,k  +PiM  l  _.P<+U  +  P/+1X+1 

Pw  —  >2  ;  >  Pe  ~  2  ’ 

.  ,  Pi,k+l  +Pi+IM  1  -  _.PiJc  +Pi+l* 

Pi  -  2  Pu~  2 


(13) 


Substitution  of  equation  (13)  into  equation  (12)  and  Use  of 
equation  (1)  lead  to 


(Pi+lk+l  -  Pi*)(Hj+\Ok  -  HjtJk+x)  +  (pt+lji  -  PiMl)(H<Ck  --ffi+lOt+l) 
Ax,Aok(Hi  +  H,+l) 


(14) 


where  Ax,-  =  x,-  +  i  —  x,  and  Act*.  =  ok  —  cs  k  +  t.  Equation  (14) 
is  the  discretization  of  the  horizontal  pressure  gradient  with 
the  finite  volume  consideration. 

[8]  Finite  difference  schemes  are  commonly  used  in  sigma 
coordinate  ocean  models.  For  the  sigma  coordinate  system, 


[9]  Discretization  of  pressure  integration  along  the  seg¬ 
ments  (equation  (11))  has  two  weaknesses:  (1)  low  accuracy 
(second-order)  and  (2)  pure  mathematical,  which  means  the 
physical  property  of  p  is  not  considered.  Since  most  sigma 
coordinate  ocean  models  are  hydrostatically  balanced, 


dp 

dz 


(16) 


one  way  to  increase  accuracy  is  to  use  Hermit  Polynomial. 

4.  Hermit  Polynomial  Integration 

[10]  Suppose  p  and  its  directional  derivative  dpldl  be 
given  at  two  vertices  of  a  segment  of  [/j,  l2 ]  of  any  finite 
volume  in  Figure  1.  Let 


=  A/=V-/,-  O7) 

varies  in  [0,  1],  The  integration  of  p  along  the  segment 
from  /]  to  l2  is  given  by, 


where 


(18) 


(!),*=+ "(!),**•  <i9> 


is  the  Hermit  polynomial  and 

=  1— 3^+2^,  4>2  =  3£,2  —  2£3, 

*i=K-2e+e,  *4=e-e, 


(20) 


are  the  four  basis  functions.  Substitution  of  equations  (19) 
and  (20)  into  (18)  leads  to 


(21) 


zi,k  =  Hi  ■  ak,  (15) 

the  horizontal  pressure  gradient  (2)  discretized  by  the 
central  difference  scheme  is 

4?\  _  Pl+iji  +PI+1XH  ~  Pit.  ~PiM  1  hk  +  o*+i\ 

Ax  Ax  2Ax,  \H,+Hm) 

,  .  >  %  *.  ■■  -  , 

(Hj+ 1  - Hj\  f Pi,k  +Pi+1/  -Pik+ 1  -Pi+u+i) 

Ax,  J  \  2Ao*  ) 

*'* )  H  * 

_  (Pl+IXH  -  Pi,k)(Hl+lak  -HjOk+ 1)  +  (pt+ix  -  Pi,k+l){HiCSk  -Hj+  l<Jk+l) 
AxiAot(Hi  +  Ht+\) 

which  is  exactly  the  same  as  equation  (14).  Thus  the 
second-order  finite  volume  scheme  is  the  same  as  the  finite 
difference  scheme  for  the  staggered  grid. 


5.  Hydrostatic  Correction 

[11]  From  the  hydrostatic  balance  (equation  (16)),  the 
pressure  integration  along  the  two  vertical  segments  of  a 
finite  volume  is  calculated  analytically  by  (taking  east  side 
in  Figure  1  as  an  example) 

[  Pedz  =  (P2  +P3)  -  {2i  ~f2g  (p2  -  p3).  (22) 

For  upper  and  lower  segments  of  the  finite  volume  (Figure  1 ), 
directional  derivative  needs  to  be  computed 


«/!/  vzi/  .  1/1/  VI/  .  .  , 

—  =  cos  a—  +  sin  a—  =  cos  a— — pgsina,  (23) 
ol  ox  oz  ox 
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[12]  The  new  correction  term  fi,*  is  in  some  sense  equiv¬ 
alent  to  the  common  practice  in  Princeton  Ocean  Model 
(POM)  [Blumberg  and  Mellor,  1987]  of  removing  a  vertical 
mean  density  profile  [ Mellor  et  al.,  1994],  though  here  is 
done  locally,  while  in  the  standard  POM  code  the  mean 
density  profile  is  usually  based  on  area  averaged  climatology 
and  requires  interpolation  of  the  mean  density  profile  to  the 
sigma  grid.  In  the  sigma  coordinate  ocean  models,  equation 
(26)  becomes 


.  '  5  ■ 

^  [(Pi+U+l  -/>■■*)  (#n-lg*  -  Hl”k+l)  +  (PMX  -PU+l)Wgt  -fli+iot+l)]  + 
-  AxiAo*(//,  +  Wj+i) 


Figure  2.  Directional  derivative  of  pressure  along  upper  or 
lower  segment. 


where  a  is  the  angle  between  the  segment  and  the  x  axis 
(Figure  2).  The  pressure  integration  along  the  upper  and 
lower  segments  is  calculated  by  (lower  segment  as  an 
example) 

•<p,-fc)  +  ^(§li -sl2)<'l“Jo)<a'2‘>' 

(24) 


Thus  integration  along  the  four  segments  of  the  finite 
volume  (Figure  1)  is  represented  by 


In  the  right-hand  side,  the  first  term  is  from  the  mean  value 
calculation,  and  the  second  term  is  the  hydrostatic  correc¬ 
tion.  Substitution  of  equation  (25)  into  equation  (10)  and 
neglect  of  high-order  terms  lead  to 


(£)„ 


[(Pn-U+1  ~  Pi,k)(zi+l,k  -  *q+l)  +  (P.+U  -Z.+l^+l)]  | 

~~  (x,+l  -  Xi)(ziJl  +  Zf+Ijt  -  Z(pt+1  -  Zf+M+l) 

(26) 


where 


£lik  = 


gF* 


6(x,+i  -  xt)  ( ziJc  +  z/+i,*  -  z.x+i  -  Zi+14+1)  ’ 
T,*  =  ^(F/f+ia*  —  //i+iCT*+i)2(pi+lt  —  Pi+i,r+i) 


/•nit 


-  -  HiOk+\)2 (Pit  -  Pi,*+i)  +  (-Hi+iut+i  -  HiOk+i )2 


The  scheme  (28)  is  evaluated  with  the  coastal  and  seamount 
topography.  The  horizontal  pressure  gradient  is  computed 
using  equations  (27)  and  (28). 


6.  Error  Estimation 

6.1.  Analytical  Coastal  Topography 

[13]  Choose  coordinates  such  that  the  y  axis  coincides  with 
the  coast,  and  the  x-axis  points  offshore.  Cross-coastal  topog¬ 
raphy  consists  of  shelf,  slope,  and  deep  layer  (Figure  3a). 
Analytical  bottom  topography  is  proposed  in  a  way  that 
shelf  and  slope  are  arcs  of  two  circles.  The  shelf  has  a 
smaller  radius  (r),  and  the  slope  has  a  larger  radius  (R).  The 
two  arcs  are  connected  such  that  the  tangent  of  the  bottom 
topography,  dh/dx,  is  continuous  at  the  shelf  break  (x  =  x0). 
This  requirement  is  met  using  the  same  maximum  expand¬ 
ing  angle  (0)  for  both  arcs  (Figure  3b).  Thus  0  represents  the 
maximum  slope  angle.  This  bottom  topography  has  three 
degrees  of  freedom:  r,  R,  and  0.  The  maximum  water  depth 
is  given  by  (Figure  3b) 


H  =  (r  +  R){  1-cosO),  (29) 


The  horizontal,  vertical  coordinates,  radii  (r,  R ),  and  water 
depth  are  nondimensionalized  by 


h*  =  ±.  (30) 


The  analytical  bottom  topography  representing  shelf,  slope, 
and  deep  layer  (Figure  3)  is  given  by 


if  (x*  <  xj) 


h*(x*)  =  < 


1  —  R*  +  tJr*1  -  (x*  -x*)2,  if  (x*  <  x*  <  xf) 

1,  if  (x*  >xf) 

(31) 


Let  the  ratio  between  the  two  radii  be  defined  by 


■  (p;+i,*+i  ~  Pa+i 


)  -  (HMak  -  HiOk)2  (pi+u  -  pa)] .  (27) 
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Figure  3.  Coastal  geometry  with  open  boundaries:  (a)  three-dimensional  view  and  (b)  cross-coastal 
view. 


Substitution  of  equations  (29)  and  (32)  into  equation  (30) 
leads  to 


the  oceanic  coordinate  system),  its  effect  on  the  sigma 
coordinate  error  (ocean)  is  neglected  in  this  study. 


(1  +  £)(!  —  cos0)’ 


x*  =  r*  sin  0  = 


R*  = 


1 


(1  -I-  £)(1  —  COS0)  ’ 
k  ■  sin0 


(1  +&)(!  —  cos0)  ’ 


6.3.  Methodology 

[is]  Error  reduction  by  the  hydrostatic  correction  term 
(Q,vt)  is  evaluated  using  the  motionless  ocean  with  an 
(33)  exponentially  stratified  density  [Chu  and  Fan,  1997], 


xf  =  (r*+R*)sin0  =  — Sm° 

1  —  cos  0 


p*  =  1  +0.005(1  -e22*),  (38) 


The  shelf  and  slope  widths  are  x*  and  (x*  -  x*), 
respectively.  Use  of  equations  (32)  and  (33)  leads  to 

x$/(xf  -  x$)  =  r*/R*  =  k  (34) 

which  indicates  that  the  nondimensional  parameter  k  is  the 
ratio  between  shelf  and  slope  widths.  Thus  the  nondimen¬ 
sional  bottom  topography  (equation  (31))  has  two  degrees 
of  freedom:  k  and  0.  Figure  4  shows  the  coastal  topography 
varying  with  the  maximum  slope  angle  0  (10°,  30°,  60°, 
90°)  and  the  shelf  slope  ratio  k  (0.1,  0.5,  1).  The  larger  the 
angle  0,  the  steeper  the  bottom  topography  is;  the  lager  the 
k,  the  shorter  the  slope  is. 


6.2.  Nondimensional  Pressure 

[14]  Let  p0  (1,025  kg  m~3)  be  the  characteristic  density. 
The  density  (p)  and  pressure  (p )  are  nondimensionalized  by 


P* 


P 

Po gH' 


(35) 


and  the  hydrostatic  balance  (equation  (16))  is  nondimensio¬ 
nalized  by 


(36) 


Integration  of  equation  (36)  from  the  ocean  surface  to  depth 
z*  leads  to 


P* 


(37) 


which  leads  to  the  fact  that  the  pressure  p*  depends  only  on 
z*  (atmospheric  pressure  effect  neglected)  and  no  horizontal 
pressure  gradient  exists.  However,  computation  in  the  sigma 
coordinate  leads  to  false  horizontal  pressure  gradient  that  is 
regarded  as  the  sigma  coordinate  error. 

[16]  The  horizontal  pressure  gradient  6p*/6x*  is  computed 
for  x*  from  0  to  12  from  the  density  field  (equation  (38)) 
with  the  analytical  bottom  topography  (equation  (31))  using 


Bottom  Topography 


where  p*  is  the  nondimensional  atmospheric  pressure.  Since  Figure  4.  Dependence  of  bottom  topography  on  k  values 
the  horizontal  atmospheric  pressure  gradient  does  not  for  different  values  of  maximum  slope  0:  (a)  10°,  (b)  30°,  (c) 
depend  on  the  ocean  bottom  topography  (or  selection  of  60°,  and  (d)  90°. 
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Figure  5.  Dependence  of  GME  on  the  maximum  slope  (0  value)  for  different  k  values:  (a)  0.1,  (b) l  0.3, 

(c)  0.6,  and  (d)  1.0.  Note  that  the  four  curves  in  each  panel  represent  the  SD  (dotted),  CD  (dot-dashed), 

CCD  (dashed),  and  HC  (solid)  schemes. 

four  different  schemes:  second-order  difference  (SD)  and  maximum  values  of  the  horizontal  pressure  gradient 
scheme  (equation  (14)),  fourth-order  compact  difference  errors. 

(CD)  scheme  [Chu  and  Fan,  1998],  sixth-order  combined  Po,. 

compact  difference  (CCD)  scheme  [Chu  and  Fan,  1998],  6-4-  Global  Performanc 

and  the  hydrostatic  correction  (HC)  scheme  (equation  (28)).  [n]  The  global  performance  of  the  four  schemes  is  - 

The  accuracy  of  these  schemes  is  evaluated  by  various  mean  uated  by  the  global  mean  (over  all  the  grid  points)  o 

r  -  /I:  v  aqok  t  •• 


Figure  6.  Dependence  of  GME  on  the  shelf  slope  ratio  (k  value)  for  different  0  value:  (a)  10°,  (b)  30°, 
(c)  60°,  and  (d)  90°.  Note  that  the  four  curves  in  each  panel  represent  the  SD  (dotted),  CD  (dot-dashed), 
CCD  (dashed),  and  HC  (solid)  schemes. 
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Vertical  Average  Pressure  Gradient 


Figure  7.  Dependence  of  VME  on  the  offshore  distance  (x*)  for  three  k  values  (0. 1 ,  0.5,  and  1 )  and  four 
0  values  (10°,  30°,  60°,  and  90°).  Note  that  each  panel  has  four  curves  representing  SD  (dotted),  CD  (dot- 
dashed),  CCD  (dashed),  and  HC  (solid)  schemes. 


absolute  values  of  the  horizontal  pressure  gradient,  denoted 
by  global  mean  error  (GME).  GME  is  computed  using  the 
four  schemes  for  various  k  (0.1  to  1.0)  and  0  (5°  to  90°). 
Figure  5  shows  the  dependence  of  GME  on  the  maximum 
slope  angle  0  with  four  different  values  of  the  shelf7slope 
ratio:  k  =  0.1,  0.3, 0.6,  1.  On  each  panel  (a  particular  value  of 
k),  four  curves  are  plotted  representing  the  four  schemes:  SD, 
CD,  CCD,  and  HC.  For  all  the  four  schemes,  GME  increases 
with  the  maximum  slope  0  monotonically.  However,  GME 
reduces  from  SD  to  CD,  from  CD  to  CCD,  and  from  CCD  to 
HC  schemes.  Taking  k  =  0.1  and  0  varies  from  5°  to  90°, 
GME  increases  from  1.04  x  10  7  to  1.58  x  10-5  using  the 
SD  scheme,  from  2.37  x  10~8  to  6.08  x  10~6  using  the  CD 
scheme,  from  1.58  x  10-9  to  5.82  x  10-6  for  the  CCD 
scheme;  and  from  1.09  x  10_u  to  1.31  x  10~8.  The  overall 
error  reduces  drastically  using  the  HC  scheme. 

[is]  Figure  6  shows  the  dependence  of  GME  on  the  shelf 
slope  ratio  (k)  with  four  different  maximum  slopes:  0  =  10°, 


30°,  60°,  and  90°.  On  each  panel  (a  particular  value  of  0), 
four  curves  are  plotted  representing  the  four  schemes:  SD, 
CD,  CCD,  and  HC.  The  same  as  Figure  5,  GME  greatly 
reduces  using  the  HC  scheme,  compared  to  using  the  three 
currently  used  schemes  (SD,  CD,  and  CCD).  For  not  very 
steep  topography  (Figures  6a-6c),  GME  reduces  from  SD 
to  CD,  from  CD  to  CCD,  and  from  CCD  to  HC  schemes. 

[19]  For  very  steep  topography  (Figure  6d,  0  =  90°),  GME 
is  comparable  using  any  of  the  existing  schemes  (SD,  CD, 
CCD),  but  2-3  orders  of  magnitude  smaller  using  the  HC 
scheme.  Such  a  high  performance  makes  the  HC  scheme 
valuable  for  coastal  modeling. 

6.5.  Cross-Coastal  Error 

[20]  The  cross-coastal  error  is  evaluated  by  the  vertical 
mean  (over  a  water  column)  of  the  absolute  values  of  the 
horizontal  pressure  gradient,  denoted  by  vertical  mean  error 
(VME).  VME  is  computed  using  the  four  schemes  for  three 
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(a)  SD:  Max:28.2  Mean:1.41  (Iff*)  (b)  CD:Max:118  Mean:3.7  (1(T7) 


Figure  8.  Contour  of  the  absolute  values  of  the  horizontal  pressure  gradient  (with  0  =  60°  and  k  =  0.2) 
calculated  using:  (a)  SD,  (b)  CD,  (c)  CCD,  and  (d)  HC  schemes. 


k  values  (0.1,  0.5,  1.0)  and  four  0  values  (10°,  30°,  60°, 
90°).  Each  panel  in  Figure  7  shows  the  dependence  of  VME 
on  die  offshore  distance  (x*)  with  four  different  schemes  for 
a  particular  combination  of  the  ( k ,  0)  values.  For  the  SD  and 
HC  schemes,  VME  usually  has  a  minimum  value  at  the 
coast  (x*  =  0)  and  increases  with  x*  from  shelf  to  slope.  For 
the  CD  and  CCD  schemes,  however,  VME  has  a  relatively 
large  value  at  the  coast,  reduces  with  x*  to  shelf  break  and 
then  increases  with  x*  in  the  slope.  This  is  caused  by 
additional  condition  such  as  cPp/dx1  =  0  is  added  at  the 
coast  (x*  =  0).  The  HC  scheme  reduces  around  105  times  of 
VME  comparing  to  the  SD  scheme,  around  104  times  of 
VME  comparing  to  the  CD  scheme,  and  around  10— 103 
times  of  VME  comparing  to  the  CCD  scheme. 

[21]  Figure  8  shows  the  contour  of  the  absolute  values  of 
the  pressure  gradient  (with  0  =  60°  and  k  =  0.2)  calculated 
using  the  four  schemes.  Large  errors  occur  near  the  bottom 
topography  for  all  the  schemes.  However,  the  error  reduces 
from  SD  to  CD,  from  CD  to  CCD,  and  from  CCD  to  HC 
schemes.  The  maximum  error  reduces  by  a  factor  of  2.39 
from  SD  to  CD  scheme,  a  factor  of  5.96  from  CD  to  CCD 
scheme,  and  a  factor  of  32.57  from  CCD  to  HC  scheme. 
GME  reduces  by  a  factor  of  3.81  from  SD  to  CD  scheme,  a 


Table  1.  Comparison  of  Maximum  and  Global  Mean  Horizontal 
Pressure  Gradient  Errors  Among  Four  Schemes  Using  the 
Analytical  Coastal  Bottom  Topography  With  0  =  60°  and  k  =  0.2 


Maximum  Error 

Global  Mean  Error 

SD 

28.2  x  1(T6 

1.41  x  10“6 

CD 

11.8  x  10~6 

0.37  x  10“6 

CCD 

1.98  x  10"6 

0.74  x  10“7 

HC 

6.08  x  10“9 

1.81  x  10“10 

factor  of  5.00  from  CD  to  CCD  scheme,  and  a  factor  of 
408.8  from  CCD  to  HC  scheme  (Table  1).  High  GME  error 
reduction  of  the  HC  scheme  is  expected  since  the  entire 
error  in  those  tests  is  due  to  the  vertical  density  gradient  (see 
equation  (38)).  Removal  of  the  hydrostatic  term  (or  hydro¬ 
static  correction)  removes  most  of  the  error. 

7.  Seamount  Test  Case 
7.1.  Known  Solution 

[22]  Accuracy  of  any  numerical  scheme  should  be  evalu¬ 
ated  by  either  analytical  or  known  solution.  For  a  realistic 
ocean  model,  the  analytical  solution  is  hard  to  find.  Consider 


Figure  9.  Seamount  geometry. 
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a  horizontally  homogeneous  and  stably  stratified  ocean  with 
realistic  topography.  Without  forcing,  initially  motionless 
ocean  will  keep  motionless  forever,  that  is  to  say,  we  have 
a  know  solution  (V  =  0).  Any  nonzero  model  velocity  can  be 
treated  as  an  error.  Ezeretal.  [2002]  evaluated  seven  different 
schemes  and  found  that  the  sixth-order  CCD  scheme  is 
accurate  but  computationally  expensive.  Thus  we  use  the 
standard  POM  seamount  test  case  to  verify  the  performance 
of  HC  scheme,  comparing  with  the  SD  and  CCD  schemes 
both  in  accuracy  and  computational  efficiency. 

7.2.  Configuration 

[23]  A  seamount  is  located  in  the  center  of  a  square 
channel  with  two  solid,  free  slip  boundaries  in  the  north 
and  south  (Figure  9).  The  bottom  topography  is  defined  by 
\Ezer  et  al.,  2002] 

H[x,y)  =H max{l  -^exp[-(x2  +/)/L2]},  (39) 


Figure  10.  Temporally  varying  MKE  (m2  s  2)  using  SD 
(dotted),  CCD  (dashed),  and  HC  (solid)  schemes. 


where  //max  is  the  maximum  depth  (4,500  m);  L  is  the 
seamount  width;  and  (l  -  A)  represents  the  depth  of  the 
seamount  tip.  In  this  study,  we  use  a  “very  steep”  case  with 
A  =  0.9  and  L-  25  km. 

[24]  Unforced  flow  over  the  seamount  in  the  presence  of 
resting,  level  isopycnals  is  an  idea  test  case  for  the  assess¬ 
ment  of  pressure  gradient  errors  in  simulating  stratified  flow 
over  topography.  The  flow  is  assumed  to  be  reentrant 
(periodic)  in  the  along  the  two  open  boundaries  of  the 
channel.  We  use  this  standard  seamount  case  of  POM  to  test 
the  performance  of  the  HC  scheme.  POM  is  the  sigma 
coordinate  model.  In  the  horizontal  directions  the  model 
uses  the  C  grid  and  the  second-order  finite  difference 
discretization  except  for  the  horizontal  pressure  gradient, 
which  the  user  has  choice  of  either  second-order  SD  or 
sixth-order  CCD  schemes.  There  are  21  sigma  levels.  For 
each  level,  the  horizontal  model  grid  includes  64  x  64  grid 
cells  with  nonuniform  size  with  finer  resolution  near  the 
seamount, 


Axi  =  Axm 


A yj  -  A yraax 


(40) 


where  Axmax  =  Aymax  =  8  km,  denoting  the  grid  size  at  the 
four  boundaries.  Temporal  discretization  is  given  by 


(At)„  =  6  s,  (A/),.„  =  30  (A/)„,  (41) 

where  (Af)e*  and  (At)in  are  external  and  internal  time  steps, 
respectively. 

[25]  POM  is  integrated  from  motionless  state  with  an 
exponentially  stratified  temperature  (°C)  field 

T(x,y,z)  =  5  +  15  exp (z/HT),  HT  =  1000  km,  (42) 

and  constant  salinity  (35  ppt)  using  three  different  schemes: 
SD,  CCD  and  HC.  During  the  integration,  a  constant 
density,  1000  kg  m-3,  has  been  subtracted  for  the  error 
reduction.  The  Laplacian  Smargorinsky  diffusion  set  up  in 
the  standard  POM  seamount  test  case  is  used  in  this  study. 


[26]  It  is  found  that  5  days  are  sufficient  for  the  model 
mean  kinetic  energy  (MKE)  per  unit  mass, 


MKE  = 


m 


pX2dxdydz 


III 


pdxdydz 


to  reach  quasi-steady  state  under  the  imposed  conditions 
(Figure  10)  for  all  the  three  schemes.  Here,  V  is  the  error 
velocity.  Thus  we  use  the  first  5  days  of  model  output  inside 
the  region  of  (-169.35  km  <  x  <  169.35  km,  -169.35  km 
<  y  <  169.35  km)  to  compare  the  HC  scheme  to  the  SD  and 
CCD  schemes.  Feasibility  of  using  the  HC  scheme  is 
twofold:  (1)  drastic  error  reduction  and  (2)  no  drastic  CPU 
time  increase.  The  following  two  quantities  are  used  to 
compare  the  errors:  the  volume-integrated  pressure  gradient 
and  the  vertically  integrated  velocity. 


7.3.  Volume-Integrated  Horizontal  Pressure  Gradient 

[27]  The  volume-integrated  horizontal  pressure  gradient 
( dp/dx )  for  a  finite  volume  ( Hdadxdy )  with  its  center  at 
(x,  y,  a),  denoted  by  VIPG,  is  used  to  represent  the  sigma 
coordinate  errors.  Figure  1 1  shows  VIPG  at  day  5  for  the 
second-order  SD,  sixth-order  CCD,  and  HC  schemes  at 
four  sigma  levels:  ka  =  5,  10,  15,  and  20  (ka  =  0  at  the 
surface,  and  ka  =  21  at  the  bottom).  VIPG  reveals  a 
dipole  pattern  at  east  and  west  of  the  seamount  at  the 
upper  half  water  column  (ka  =  5,  10),  and  a  quadruple 
dipole  in  the  lower  layer.  VIPG  increases  with  depth  from 
the  surface  (maximum  3.41  x  104  N)  to  the  bottom 
(maximum  3.98  x  105  N)  using  the  SD  scheme).  It  reduces 
greatly  using  both  CCD  and  HC  schemes.  Taking  ka  =  10  as 
an  example,  the  maximum  VIPG  is  2.01  x  105  N  using  the 
SD  scheme,  3.42  x  103  N  using  the  sixth-order  CCD 
scheme,  and  3.05  x  103  N  using  the  HC  scheme.  This 
indicates  that  the  HC  scheme  has  a  comparable  accuracy  as 
the  sixth-order  CCD  scheme  and  reduces  the  errors  by  a 
factor  of  70  comparing  to  the  SD  scheme. 

7.4.  Vertically  Integrated  Velocity 

[28]  Owing  to  a  very  large  number  of  calculations 
performed,  we  discuss  the  results  exclusively  in  terms  of 
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Figure  11.  VIPG  (unit:  N)  on  day  5  at  four  sigma  levels:  ka  =  5,  10,  15,  and  20  (ka  =  0  at  the  surface, 
and  ka  =  21  at  the  bottom)  using  SD,  CCD,  and  HC  schemes. 


the  vertically  integrated  velocity  (Vi„t)  generated  by  the 
pressure  gradient  errors.  Figure  12  shows  the  time  evolu¬ 
tion  of  the  vertically  integrated  error  velocity  for  day  1, 
day  2,  day  4,  and  day  5  using  the  SD,  CCD,  and  HC 
schemes. 

[29]  The  vertically  integrated  velocity  (Vint)  has  a  large- 
scale  eight-lobe  pattern  centered  on  the  seamount.  This 
symmetric  structure  can  be  found  in  all  the  fields.  However, 
its  maximum  values  are  around  0.05-0.1  m2  s-1  using  the 
SD  scheme,  0.001-0.002  m2  s-1  using  the  CCD  and  HC 


schemes.  After  5  days  of  integration,  the  model  generates 
spurious  vertically  integrated  currents  of  0(0.09  m2  s-1) 
using  the  SD  scheme  and  of  0(0.0015  m2  s-1)  using  the 
CCD  and  HC  schemes.  Thus  the  HC  scheme  has  the  same 
accuracy  as  the  sixth-order  CCD  scheme. 

7.5.  Computational  Efficiency 

[30]  The  numerical  integration  of  POM  using  the  three 
schemes  was  performed  using  SGI-Origin  200  machine. 
The  CPU  time  was  recorded  at  the  end  of  day  5  integration. 
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Figure  12.  Vertically  integrated  velocity  vectors  (unit:  m2  s  ')  at  day  1,  day  2,  day  4,  and  day  5  using 
SD,  CCD,  and  HC  schemes. 


It  increases  36%  from  the  SD  to  CCD,  and  only  3%  from 
the  SD  to  HC  scheme  (Table  2).  No  evident  CPU  increase 
from  the  SD  to  HC  scheme  but  the  drastic  error  reduction 
make  the  HC  scheme  a  good  choice  for  the  sigma  coordi¬ 
nate  ocean  model  to  reduce  the  horizontal  pressure  gradient 
error. 

8.  Comparison  Between  the  Two  Test  Cases 

[31]  Comparing  to  the  analytical  coastal  topography  test, 
the  seamount  test  case  yields  smaller  error  reduction.  The 
large  difference  in  the  error  reduction  is  caused  by  the 
existence  of  the  two  kinds  of  the  sigma  error  [Mellor  et  al., 
1998].  The  two-dimensional  coastal  topography  test  only 


includes  the  sigma  error  of  the  first  kind,  which  often 
decrease?  with  time.  While  the  three-dimensional  seamount 
test  case  also  includes  the  sigma  error  of  the  second  kind 
that  depends  on  the  curvature  of  the  topography  and  does 
not  reduce  with  time.  Thus  high-order  interpolation 


Table  2.  Comparison  of  CPU  Time  (Minute)  at  the  End  of  5  Day 
Runs  of  the  POM  Seamount  Test  Case  Using  the  SD,  CCD,  and 
HC  Schemes 


Scheme 

SD 

CCD 

HD 

CPU  Time 

171.51 

233.33 

176.92 

Ratio 

1 

1.36 

1.03 

37-12 
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schemes  are  especially  useful  for  the  second  kind  of  the 
sigma  error. 

9.  Conclusions 

[32]  1.  The  sigma  coordinate  pressure  gradient  error 
depends  on  the  choice  of  difference  schemes.  By  choosing 
an  optimal  scheme,  we  may  reduce  the  error  significantly 
without  increasing  the  horizontal  resolution.  An  optimal 
scheme  should  satisfy  the  three  requirements:  conservation 
(especially  for  pressure  gradient  calculation),  high  accuracy 
(especially  near  the  steep  topography),  and  computational 
efficiency. 

[33]  2.  The  hydrostatic  correction  scheme  has  been  de¬ 
veloped  in  this  study  using  the  finite  volume  discretization. 
The  major  features  of  this  new  scheme  are  conservation  (for 
pressure  gradient  calculation),  high  accuracy,  and  computa¬ 
tional  efficiency. 

[34]  3.  Computation  of  horizontal  pressure  gradient  for 
analytical  coastal  topography  with  two  varying  parameters 
(maximum  slope  and  shelf  slope  ratio)  is  used  to  verify  the 
performance  of  the  hydrostatic  correction  scheme  against 
the  second-order,  fourth-order  combined,  and  sixth-order 
combined  compact  schemes.  For  all  the  four  schemes,  the 
global  mean  error  increases  with  the  maximum  slope  and 
decreases  with  the  shelf  slope  ratio  monotonically.  Howev¬ 
er,  it  reduces  by  a  factor  of  3.81  from  the  second-order  to 
fourth-order  compact  schemes,  a  factor  of  5.00  from  the 
fourth-order  compact  to  sixth-order  combined  compact 
schemes,  and  a  factor  of  408.8  from  the  sixth-order  com¬ 
bined  compact  to  hydrostatic  correction  schemes.  The 
maximum  error  reduces  by  a  factor  of  2.39  from  the 
second-order  to  fourth-order  compact  schemes,  a  factor  of 
5.96  from  the  fourth-order  compact  to  sixth-order  combined 
compact  schemes,  and  a  factor  of  32.57  from  the  sixth-order 
combined  compact  to  hydrostatic  correction  schemes. 

[35]  4.  For  very  steep  topography,  the  global  mean 
pressure  gradient  error  is  comparable  using  any  of  the 
existing  schemes  (second-order,  fourth-order  compact, 
sixth-order  combined  compact),  but  2-3  orders  of  magni¬ 
tude  smaller  using  the  hydrostatic  correction  scheme.  Such 
high  performance  of  the  hydrostatic  correction  scheme  is 
due  to  its  conservation  characteristics  for  calculating  the 
horizontal  pressure  gradient. 

[36]  5.  The  standard  POM  seamount  test  case  is  used  to 
evaluate  the  performance  of  the  hydrostatic  correction 
scheme  against  the  second-order,  and  sixth-order  combined 
compact  schemes.  The  hydrostatic  correction  scheme  has 
comparable  accuracy  as  the  sixth-order  combined  compact 
scheme  and  reduces  the  errors  by  a  factor  of  70  comparing 
to  the  second-order  scheme.  However,  its  computational 
efficiency  is  comparable  to  the  second-order  scheme. 


[37]  6.  As  model  resolution  increases,  the  cost  for  given 
accuracy  will  eventually  favor  the  high-order  methods. 
While  the  HC  scheme  looks  the  best  overall  here,  more 
stringent  accuracy  requirements  could  be  easily  satisfied 
using  the  HC  scheme. 

[38]  Acknowledgments.  The  Office  of  Naval  Research,  the  Naval 
Oceanographic  Office,  and  the  Naval  Postgraduate  School  supported  this 
study. 
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